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ABSTRACT 



Context. Observations of coUimated outflows in young stellar objects indicate that several features of the jets can be understood by 
adopting the picture of a two-component outflow, wherein a central stellar component around the jet axis is surrounded by an extended 
disk-wind. The precise contribution of each component may depend on the intrinsic physical properties of the YSO-disk system as 
(~| , well as its evolutionary stage. 

Aims. In this context, the present article starts a systematic investigation of two-component jet models via time-dependent simulations 
I ' of two prototypical and complementary analytical solutions, each closely related to the properties of stellar-outflows and disk-winds. 

O ' These models describe a meridionally and a radially self-similar exact solution of the steady-state, ideal hydromagnetic equations, 

respectively. 

C/3 , Methods. By using the PLUTO code to carry out the simulations, the study focuses on the topological stability of each of the two 

^ , ■ analytical solutions, which are successfully extended to all space by removing their singularities. In addition, their behavior and ro- 

bustness over several physical and numerical modifications is extensively examined. Therefore, this work serves as the starting point 
for the analysis of the two-component jet simulations. 
^ . Results. It is found that radially self-similar solutions (disk-winds) always reach a final steady-state while maintaining all their well- 

' defined properties. The different ways to replace the singular part of the solution around the symmetry axis, being a first approximation 

^— ■») ' towards a two-component outflow, lead to the appearance of a shock at the super-fast domain corresponding to the fast magnetosonic 

, separatrix surface. These conclusions hold true independently of the numerical modifications and/or evolutionary constraints that 

the models have been undergone, such as starting with a sub-modified-fast initial solution or different types of heating/cooling as- 
sumptions. Furthermore, the final outcome of the simulations remains close enough to the initial analytical configurations showing 
thus, their topological stability. Conversely, the asymptotic configuration and the stability of meridionally self-similar models (stellar- 
winds) is related to the heating processes at the base of the wind. If the heating is modified by assuming a polytropic relation between 
, density and pressure, a turbulent evolution is found. On the other hand, adiabatic conditions lead to the replacement of the outflow by 

■ an almost static atmosphere. 



Key words. ISM/Stars; jets and outflows - MHD - Stars: pre-main sequence, formation 
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' 1 . Introduction on the different driving mechanisms proposed, it is not yet clear 

which is the dominant plasma launching mechanism in YSO jets. 

Observations made over the last two decades have shown that ^j^^ ^ ^ protostellar object basically contains two 

one class of the widespread astrophysical phenomenon of colh- dynamical constituents, a central protostar and its suiTounding 

mated plasma outflows (]ets) is being launched from the vicin- ^^cretion disk. Consequently, in Bogovalov & Tsinganos OOOl) 

ity of most young stellar objects (YSOs) (Burrows et al. [19%^. ^ ^^^^ -^^^ observed from T Tauri stars most likely 

These supersonic mass ot^ows are found to becorrelated with ^^^^-^^ ^^^^^ components: (i) an inner pressure 

accretion (Cabrit et al. [IMl Hartigan et al. [1995), to have nar- ^^-^^^ ^-^^^ ^^^-^^ non-collimated if the star is an inefficient 

row opening angles (Ray et al. fl996^ and to propagate for sev- j^agnetic rotator and (ii) an outer magneto-centrifugally driven 

era! orders of magnitude of spatial distances ranging from the jisk-wind which provides most of the high mass loss rate ob- 

AU to the pc scales (Dougados et al. 2000 ; Hartigan et al. 2004f ^^^^^^ relatively faster rotating magnetized disk produces 

A central role m the launching, acceleration and colliination of ^j^^ self-collimated wind which then forces all enclosed outflow 

-"^/AJ^x^^^^ believed to be played by niagnetohydrody- ^^^^ ^j^^ ^^^^^^^ ^^^^^^ collimated as well. This conclusion 

namic (MHD) eff^ects, which can also successfully remove the confirmed by self-consistent simulations of the MHD equa- 

excessive angular momentum, allowing in this way the YSO to ^-^^^ ^^^.^ ^^^^^^^ ^^^^^-^^ ^1 ^ ^j^^^ 

accrete and enter the inain sequence. Nevertheless, although re- ySO jets observed in association with T Tauri stars, in 

cent high angular resolution observations put several constraints ^^^^-^^ p^^^^^^.^ ^^-^^^ ^^^U^ ^^^^^^ magneto- 

centrifugally launched extended "warm" disk-wind, a third com- 

Send offprint requests to: Titos Matsakos, ponent may be driven by magnetic processes at the magneto- 
e-mail: matsakos(aph. unite, it sphere/disk interaction, i.e., a sporadically ejected X-type wind. 
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In addition, in Ferreira et al. (|2000| l a non steady "two-flow" 
scenario was also suggested, regarding a reconnection X- and a 
disk-wind. Nevertheless, the existence of such sporadic compo- 
nents is not supported by observational data as being one of the 
major contributors to the steady characteristics of jets, but rather 
could explain the observed variability in jet emission. On the 
other hand and within the same framework, recent observations 
(Edwards et al. 2006; Kwan et al. 2007 ) show that both, disk 
and stellar winds are mainly present in T Tauri stars, with the 
dominant component being determined by the intrinsic physical 
properties of the particular YSO. 

From the analytical studies point of view, the complexity of 
the launching and coUimation mechanisms of jets have forced 
researchers over the past several years to treat these two com- 
ponents separately. The only available analytical MHD models 
for jets are those characterized by the symmetries of radial and 
meridional self-similarity (Vlahakis & Tsinganos |1998l l. In the 
former case, the solution is invariant as we look at a constant po- 
lar angle and in the latter, as we look at a constant spherical ra- 
dius. The computational consequence of the respective symme- 
try is that by employing the separable spherical coordinates (r, 
ff), the set of coupled MHD equations reduces to a set of ordinary 
differential equations in 6, or, in r, respectively. The last remain- 
ing difficulty is to select solutions which are causally discon- 
nected from the source of the outflow, i.e., those crossing the fast 
magnetosonic separatrix. In this way, one may construct either 
radially self-similar solutions closely related to the properties 
of magneto-centrifugally driven disk- winds (Blandford & Payne 
LL982; Contopoulos & Lovelace 1994; Ferreira 1997; Vlahakis et 
al. |2000 hereafter VTSTOO), or meridionally self-similar ones 
to address pressure driven stellar outflows (Sauty & Tsinganos 
[19941 Trussoni et al. [19971 Sauty et al. [20021 hereafter STT02). 

Since each self-similar symmetry corresponds to a particular 
component, we adopt the following initials: ADO (Analytical 
Disk Outflow) and ASO (Analytical SteUar Outflow) to refer to 
radially and meridionally self-similar solutions, respectively. 

Apart from the geometry, an intrinsic distinction between 
these two classes concerns the treatment of the energy equation. 
Nevertheless, the symmetry difference makes them complemen- 
tary to each other, since the ADO solution becomes singular at 
the axis, whereas the ASO is by definition the proper one for 
modeling the area close to it. In addition, the properties of the 
launching region of the disk-wind, i.e., at large polar angles, are 
described more naturally by the ADO model. For a recent review 
on the analytical work on MHD outflows the reader is referred 
to Tsinganos (I2007I I. 

On the other hand, the increase of computational power 
along with the development of sophisticated numerical codes 
have allowed to study the time evolution of the MHD equations, 
giving us new perspectives of the physics involved. Jet launching 
and collimation have been mainly investigated with the follow- 
ing two methodologies: (i) by treating the disk as a boundary 
(Krasnopolsky et al. ,1999. Ouyed et al. ,2003. Fendt .2006) and 
(ii) by including the disk inside the computational box (Casse & 
Keppens'2004', Meliani et al. [20061 Zanni et al. [2007l l. The for- 
mer case allows a wider range of physical processes and mecha- 
nisms to be studied, whereas the latter, has the advantage of the 
jet evolution being consistent with that of the disk. However, 
with the exception of Gracia et al. (120061 hereafter GVT06), 
most numerical studies did not take advantage of the availability 
of the well studied analytical solutions which also allow a para- 
metric study and therefore a better physical understanding of the 
problem of jet launching and collimation. 



The present work is the first attempt to numerically construct 
and study a two-component jet, by using as starting point the two 
well studied classes of analytical self-similar solutions. Towards 
this goal, in this paper we first address the question of the topo- 
logical stability of each one of these two classes separately, be- 
fore we combine them in the following paper 

Concerning the ADO model, we shall use the VTSTOO an- 
alytical solution, completing and considerably extending the 
GVT06 analysis. Therein, they found that the disk-wind model 
may attain a new steady-state configuration close enough to the 
initial analytical one, provided the assumption of some appropri- 
ate approximations around the axis. We further present the first 
numerical studies of ASO (meridionally self-similar winds), re- 
ferring to the solutions of STT02, that are essential to model the 
region around the axis, just where the ADO model fails. 

Once the physical properties and stability of these two 
classes of solutions is clarified, our final aim will be to effec- 
tively build up a model that consistently merges the ASO and 
ADO solutions. Such simulations will be presented in a future 
work where we will study the launching and propagation of a 
collimated stellar wind around the system axis surrounded by a 
disk wind. 

Finally, a word on the term topological stability used in this 
paper. Classical stability theory addresses the question whether a 
given equilibrium configuration evolves away from (=unstable) 
or back to (=stable) the initial equilibrium when perturbed. In 
the present context, topological (or structural) stability refers to 
the question whether a given configuration preserves its topolog- 
ical properties when subject to various perturbations. Needless 
to say, that topologically unstable configurations may well be 
stable from the classical point of view and vice versa. 

The paper is structured as follows: In section ^ the for- 
malism of the ADO and ASO solutions is briefly presented. In 
section ^ their implementation is explained and the numeri- 
cal models to be investigated are presented. Section ^reports 
the results obtained by carrying out the respective simulations. 
Finally in section ^we discuss our results in the framework of 
the future matching and we report the conclusions of this work. 



2. MHD equations and the self-similar solutions 

The ideal MHD equations are: 



^ + v-(py) = o, 

ot 

dV 1 1 

— + (V ■ V)V + -B X (V X B) + -VP = -VO , 

ot p p 

— + V-VP + Ff V ■ y = A , 

dt 

dB 

— -Vx(yxB) = 0, 
ot 



(1) 

(2) 
(3) 
(4) 



where p, P, V, B denote the density, pressure, velocity and mag- 
netic field over V47r, respectively. O = -QMjR is the grav- 
itational potential of the central object {Q is the gravitational 
constant) with mass M, A represents the volumetric energy 
gain/loss terms (A = [F - l}pQ, with Q the energy source terms 
per unit mass), and F is the ratio of the specific heats. 

For the sake of clarity, we adopt the following notation: the 
subscripts r and 9 of the physical variables, will be used to re- 
fer to the ADO and ASO solutions respectively, whereas the 
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cylindrical radial direction is denoted with the symbol m and 
the spherical radial direction with R. Finally, the index p corre- 
sponds to the poloidal components of the physical variables. 

By assuming steady-state and axisymmetry, several con- 
served quantities exist along the fieldlines (Tsinganos il982i) . By 
introducing the magnetic flux function A - ( j^) f Bp • c/S to la- 
bel the iso-surfaces that enclose constant poloidal magnetic flux, 
then these integrals take the following simple formulation: 



(5) 



L(A) 



pVp 




" Bp 












m \ 


P 




BA 





(7) 



where is the mass-to-magnetic-flux ratio, Q. the field angular 
velocity, and L the total specific angular momentum. The ratio 
y[LjTi defines the Alfvenic lever arm at each fieldline, where 
the poloidal flow speed is equal to the poloidal Alfvenic one. In 
the adiabatic-isentropic case where A = 0, there exist two more 
integrals, the total energy flux density to mass flux density E and 
the specific entropy Q, which are given by; 



y2 

E{A) = — -H - 

^ 2 r 



G(A) = 4 ■ 



P 



(8) 



(9) 



Hence, one obtains the following expressions for the physical 
variables: 



A--3/2 



1 



Pr = PrX 



■-2 



2y 



, ,4 Ml sin 
V - -V a ' - 



(6) V^,r = V,.,Aa; 



1/4, 



G2 cos(iA,. + 0) 
G,(l-M2)' 



(cos + sin i//,-z) 



B - -R (7^'^^"' 



sin 9 



1 -G? 
G,(l-M2)- 



(cos if/r-uj + sin l//rZ) , 



(13) 
(14) 

(15) 
(16) 
(17) 

(18) 



The starred quantities are related to their characteristic values 
at the Alfven radius m, along the reference fieldline a = I. 
Moreover, they are interconnected with the following relations: 



Br, 

^Ip7* 



Pr, = 



0M 



(19) 



where the constants A measures the strength of rotation and TC 
the gravitational potential. Finally // is proportional to the gas 
entropy. 



2.1. ADO - The radially self-similar model 

We employ the radially self-similar solution which is described 
in VTSTOO and crosses successfully all three critical surfaces. 
We note that a polytropic relation between the density and the 
pressure is assumed, i.e. P - 2(A)p^, with y being the effective 
polytropic index. Equivalently, the source term in Eq. (|3]l has the 
special form 



A = (F - r)P(V ■ V) , 

transforming the energy equation ([3]) to 

— + V ■VP + yPV ■ y = . 



(10) 



(11) 



The latter can be interpreted as describing the adiabatic evolution 
of a gas with ratio of specific heats y, whose entropy Pjp'^ is 
conserved. 

The solution is provided by the values of the key functions 
Mr{0), Gr{9) and iAj-(^?), which are the Alfvenic Mach number, 
the cylindrical distance in units of the corresponding Alfvenic 
lever arm and the angle between a particular poloidal fieldline 
and the cylindrical radial direction, respectively. Then, the field- 
lines can be labeled byQ 



Br*tjji 



where o-,- 



\Gl 



(12) 



' Note that x is a model parameter governing the scaling of the mag- 
netic field and is related to ^ = 2(x - 3/4) which is a local measure of 
the disk ejection efliiciency in the disk model of Ferreira I il997t . 



2.2. ASO - The meridionally self-similar model 

We employ a meridionally self-similar solution which corre- 
sponds to the case of a spherically symmetric thermal pressure 
(case K = 0, second sub-table of Table 1 in STT02, where k 
represents the deviations from such a pressure symmetry). This 
solution is derived without the assumption of a polytropic rela- 
tion, with the respective energy source term being consistently 
derived a posteriori. In this case the key functions are Gg(R), 
Mg(R), Hg{R) and Fg(R). The former two have the same interpre- 
tation as in the ADO model, while the latter ones are the pres- 
sure and the expansion factor, respectively. The magnetic flux A 
is then given by 



, Be*Ri 



where ag - 



RlG^ 



(20) 



Moreover, the physical variables are given from the follow- 
ing expressions: 



pe - pg*—^{l + Sag) , 



Pe = i'e.nfl , 

Ml sine cose 



G^ VI + Sag 



('-?)■ 



IVlg 



' Gi VI + Sag 



H- dag ^ 



cos^ + sin^ 



Fg\ 



(21) 
(22) 
(23) 

(24) 
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(-1 

Vd,:0 = VeU a, 



■Mi 



1 



Gg(l-Ml) VI + dae ' 
sin cos / Fg\ 

^(cos^e + sin^e:^) , 



Bs-e - —BetA'aJ^- 



1-Gl 



(25) 
(26) 
(27) 
(28) 



Here, 6 describes deviations from a spherically-symmetric den- 
sity whereas A' the strength of the magnetic torque at the Alfven 
radius The starred quantities are the reference values at 
and are related as follows: 



Bg 



1 7 
-Bl 



29M 



(29) 



where v represents the strength of the gravitational potential. 
Finally, the expression for the energy source term is 



A 



PeVR;g Vl 



1 -I- dag 2Rt 



,dU„ 



dM; 



2\ 



Mi — - + rng 

^ dR dR ) 



(30) 



2.3. Physical aspects and differences of the two solutions 

In order to have a better understanding of these two classes of 
self-similar solutions, we pause here to point out a few aspects 
concerning the main physical mechanisms involved in each case 
and their major intrinsic differences. 

The ADO (radially self-similar) solution corresponds to a 
magneto-centrifugally driven outflow with the "bead on a ro- 
tating wire" analogy. Fig. 7 of VTSTOO displays the different 
terms of the conserved total energy as plotted along a particular 
fieldline as a function of z. It is evident, that close to the base of 
the outflow, the electromagnetic energy dominates, whereas, as 
the flow is being accelerated, the poloidal kinetic one becomes 
eventually the main component of the total energy. The minor 
role that the enthalpy seems to play is being investigated in ^ 
Finally the reader is also referred to Fig. 5 of VTSTOO where the 
components of the outflow speeds are plotted. 

On the contrary, the ASO (meridionally self-similar) solution 
adopted is a pressure driven outflow. Fig. 9 of STT02 presents 
the forces acting along and across the streamlines for a solution 
very similar to the one employed here. It can be clearly seen 
that the pressure gradient is the dominant force for the acceler- 
ation of the flow while the collimation mechanisms are due to 
the hoop stresses. Another important feature of the ASO model 
is that apart from the polar fieldlines leaving the stellar sur- 
face and closing to infinity, there are also those which cross the 
equator being asymptotically parallel to the axisymmetry axis. 
Moreover, a "dead-zone" also exists and is defined by the re- 
gion of the fieldlines with both their footpoints rooted on the 
star (Fig. 10 of Sauty & Tsinganos |1994| ). Furthermore, we note 
that meridionally self-similar models are classified by an ener- 
getic criterion which characterizes the asymptotic shape of the 
streamlines. In our case, the employed ASO corresponds to a 
cylindrical, magnetically collimated jet. 

A common feature concerning both classes, is the fact that 
the poloidal critical surfaces do not coincide with the surfaces 
where the steady-state MHD equations change character from 



elliptic to hyperbolic and vice versa. This is because of the con- 
straint put on the propagation of the MHD waves by the axisym- 
metry and the self-similarity assumption (Tsinganos et al. 19961. 

As far as the two main intrinsic differences of the two models 
are concerned, we note the following. First, due to the assump- 
tions of self-similarity, the ADO model has its MHD critical sur- 
faces given for 6cr = const, hence being of conical shape. On the 
other hand, they are spherical in the case of the ASO, as an out- 
come of the radial dependence of its key functions. 

The second difference concerns the energy equation. In the 
ADO solution, by assuming a polytropic relationship between 
density and pressure, the total energy-to-mass-flux-ratio is con- 
served. However, in the ASO, the momentum equation pro- 
vides enough relations to close the system and the total energy- 
to-mass-flux-ratio is not used. In both cases though, the heat- 
ing/cooling mechanisms necessary to maintain the outflow can 
be calculated a posteriori (see Eq. fTOl and 1301 ). 

3. The numerical models 

We will mainly focus on four topics sorted by increasing impor- 
tance: 

1 . to complete the GVT06 work by imposing the correct num- 
ber of boundary conditions according to the number of waves 
propagating downstream; 

2. to further extend the GVT06 results by including the equa- 
tor inside the computational domain and by investigating the 
effects of the singularity substitution, the resolution and the 
choice of the minimum vertical distance of the lower bound- 
ary of the computational box, z„„„ ; 

3. to carry out and present the first time-dependent simulations 
of an ASO solution; 

4. to study how the energy input/output modifications influence 
the features, stability and robustness of each model. 

The analytical solutions provide the key functions, akeady 
discussed, along 6 and R for the ADO and ASO models, respec- 
tively. Then, by properly interpolating in a cylindrical or a spher- 
ical grid the physical values are initialized with the help of Eqs. 
(fT2]i-(fT8Tl and (E0b-(|28]). More details on the treatment of the axis 
concerning the ADO solution are given in section ^ 

Eqs. ([T]i-@ are solved numerically using the MHD module 
provided by the PLUTO code0 (Mignone et al. 2007 ). PLUTO is 
a modular Godunov-type code particularly oriented towards the 
treatment of astrophysical flows in the presence of discontinu- 
ities. For the present case, second order accuracy is achieved us- 
ing a Runge-Kutta scheme (for temporal integration) and piece- 
wise linear reconstruction (in space). Although all the compu- 
tations were carried out with the simple (and computationally 
efficient) Lax-Friedrichs solver, we point out that no significant 
differences were found by switching to more complex Riemann 
solvers available in the code. 

Table [1] lists the numerical models constructed in order to 
study the previously mentioned aspects. Note that the ADO 
model is being investigated more extensively due to the singu- 
larity appearing at small polar angles. We define the reference 
lengths -Wt and R,, of the ADO and the ASO models respec- 
tively, to be unity. In addition, the reference velocities are nor- 
malized by setting V,-, - 1 and Vg, - V2"7C/v in order for both 
solutions to have the same gravitational potential. Time will be 

expressed in units of f, = In ^j■uJl|QM - 2n '<JmlVrt/'K - n, 



Publicly available at http://phitocode.to.astro.it 
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Table 1. List of the numerical models. The first 14 lines refer to Analytical Disk Outflows (ADO) solutions, the last 7 lines to the 
Analytical Stellar Outflow (ASO) solutions. 



Name 


Model 


Geometry 


Grid [m x z] or [R x 9] 


Resolution 


Total time 


Description 


WR 


ADO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


6.0 


GVT06 (overspecified b.c.) (Fig.|2|l 


2DB 


ADO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


6.0 


Correct number of b.c. (Fig.|2l[8]l 


3DB 


ADO 


Cylindrical 


[0, 50] X [6, 100] 


256 X 480 


6.0 


Correct number of b.c, higher resolution, (Fig.[Tll 


4DH 


ADO 


Cylindrical 


[0, 50] X [6, 56] 


400 X 400 


2.0 


Shock study 


SDH 


ADO 


Cylindrical 


[0,50] X [6,56] 


800 X 800 


2.0 


Shock study (Fig.Elllll 


6DH 


ADO 


Cylindrical 


[0, 50] X [6, 56] 


1200 X 1200 


2.0 


Shock study (Fig.|3j 


IDS 


ADO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


6.0 


Different singularity smoothening (Fig.|6l[8]l 


SDS 


ADO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


6.0 


Sub modified fast lower boundary (Fig.|6l|7][8j 


9DZ 


ADO 


Cylindrical 


[0, 50] X [3, 100] 


128 X 248 


15.0 


Lower z„„„ (Fig.|9j 


lODZ 


ADO 


Cylindrical 


[0,50] X [12, 100] 


128 X 224 


4.0 


Higher z„„„ (Fig.lgJ 


IIDE 


ADO 


Spherical 


[10,90]x[0,;r/2-6] 


408 X 200 


2.0 


Extension up to the equator fFig.llOt 


ilDG 


ADO 


Spherical 


[10,90] X [0,;r/2- 6] 


408 X 200 


10.0 


Adiabatic, F = 5/3 (Fig.fTOt 


13DI 


ADO 


Spherical 


[10,90] X [0,n/2-E] 


408 X 200 


2.0 


Isothermal (Fig.llOil 


14DT 


ADO 


Cylindrical 


[0,200] X [6, 100] 


512 X 240 


220.0 


Long term simulation (Fig. 11 111 


ISR 


ASO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


4.0 


Super-Alfvenic domain in cylindrical fFig.ll2lll4t 


2SL 


ASO 


Spherical 


[1, 1000] X [0,;r/2- 6] 


200 X 200 


40.0 


Log. grid, super- Alfvenic domain 


3SL 


ASO 


Spherical 


[0.7,7] X [0,7r/2- e] 


200 X 200 


4.0 


Log. grid, sub- Alfvenic included (Fig.ll3|[T5t 


4SP 


ASO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


600.0 


Use of polytropic index y = 1.05 (Fig.ll4t 


5SG 


ASO 


Cylindrical 


[0, 50] X [6, 100] 


128 X 240 


600.0 


Adiabatic, F = 5/3 (Fig.[T4j 


6SL 


ASO 


Spherical 


[0.7,7]x [0,;r/2-6] 


200 X 200 


80.0 


Log. grid, trans-Alfvenic, polytropic (Fig. 11 6b 


ISL 


ASO 


Spherical 


[0.7,7]x [0,7r/2-6] 


200 X 200 


80.0 


Log. grid, trans-Alfvenic, adiabatic (Fig.ll6t 



i.e. using the Keplerian period at distance or on the equa- 
torial plane. The final time of the simulations is obviously cho- 
sen to be greater than the one needed for a steady-state to be 
reached, if of course there exists one. Notice that all the fol- 
lowing figures, apart from Figs. |6] and [16] correspond to the fi- 
nal state reached by the simulation at the final time indicated 
in Table [T] In order to prove the stability of the steady state 
we have included a long term simulation to make the argument 
solid. Finally, concerning the choice of our computational box, 
in many models we follow the guidelines of GVT06. 

3.1. The initial ADO model 

The model parameters of this solution were chosen as x = 0.75 
and 7 = 1.05, while the solution parameters are given tohe A - 
11. 70, fi = 2.99, = 2.00, in accordance to VTSTOO0 

By definition, radially self-similar models fail to provide 
physically accepted solutions close to the axis, due to the lo- 
cal diverging behavior of the physical variables. The solution of 
VTSTOO is terminated at 6* < 6I„„„ = 0.025(rad) - 1.5°, after 
having crossed the modified fast critical surface, also known as 
the fast magnetosonic separatrix surface (FMSS) (Tsinganos et 
al. |1996l l. However, in GVT06 the rotational axis was included 
self-consistently inside the computational box by properly ini- 
tializing the physical variables at these small polar angles. This 
was achieved by assuming that the key functions G{d), M(ff), and 
i//{d) are all even, e.g. G{6) - G{-ff), and then by accordingly in- 
terpolating. The time evolution was performed with the numer- 
ical code NIRVANA. We have implemented this setup (model 
IDR) in PLUTO, although with a slightly different extrapolation 
scheme. Our proper smoothening of the flow quantities near the 



^ The solution adopted with the model parameter x = 0.75 corre- 
sponds to a zero ejection index f according to Ferreira (T997i. However, 
as it is evident from Figs. 5 and 6 of VTSTOO the solution with 
X = 0.7575, i.e. f = 0.0025, is almost identical to the one with x = 0.75 
for z > 0. 1 . Therefore, we argue that the ADO solution employed here 
should not contradict the theoretical arguments presented in Ferreira 



axis follows the same guidelines, i.e. linearly extrapolating the 
key functions and then initializing the physical variables. This 
is discussed in detail later on, where we apply and investigate 
the effect of other extrapolation schemes as well. Finally, notice 
that such a modification of the inner part of the wind, mimics the 
presence of an eff'ective stellar outflow. 

Obviously, such assumptions would not retain a divergence- 
free magnetic field in the central region, if Eq. ( [TtT i is used for its 
initialization. So, instead, we initialize the poloidal component 
of the magnetic field by taking advantage of the flux function 
(fT2] i. Finally, in order to ensure consistency with the ideal MHD 
assumption, the radial velocity component is derived from the 
relation: 

V^^B^^. (31) 
where is taken from Eq. (fTSl l. 
3.2. The initial ASO model 

The adopted analytical solution has the model parameter /c = 0, 
and corresponds to the following values of the solution parame- 
ters: 6 = 0.01, ^' = 3, V = 24.1, in accordance to STT02. 

In this case, the values of the key functions are available 
for Q.6R* < R < IQ'^R*. However, note that this model varies 
significantly over all its radial scales and hence it is numeri- 
cally complicated to resolve all regions with adequate accuracy. 
For this reason, a number of simulations consider the super- 
Alfvenic domain only. On the other hand, exploiting the fact 
that PLUTO can integrate in time over a uniform logarithmic 
grid, the full range of the solution is also studied. Furthermore, 
the energy source term, given by expression ( [30l l. is being taken 
into account during the numerical evolution, unless otherwise in- 
dicated. Notice however, that the energy source term is not free 
to evolve in time, since it is provided by the key functions of 
the analytical solution and hence it is kept fixed throughout the 
simulations. 
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3.3. Boundaries 

The number of boundary conditions imposed at a physical 
boundary has to match the number of characteristic waves car- 
rying information from the boundary towards the inside of the 
computational box (Bogovalov 1997). This defines the physical 
boundary conditions (Thompson 19871 |l990t Mignone l2005l l as 
opposed to the information directed outside, which can be en- 
tirely determined by the solution inside the computational box. 
Still, since our numerical scheme requires the knowledge of all 
(eight) flow variables in the ghost zones, additional numerical 
boundary conditions are prescribed using suitable one-sided ex- 
trapolation formula. Furthermore, although the physical quanti- 
ties evolved in the code are 8, i.e., p, P, B, V, axisymmetry along 
with the V ■ B = condition reduces the number of variables to 
7. 

Note that in GVT06, all the quantities were kept fixed at 
the lower z boundary postponing the examination of this over- 
specification issue for a future study. 

In this framework, when cylindrical coordinates are adopted, 
we divide the lower boundary in 4 regions, i.e. a) V, > V/oj;;;, b) 

Vfast,z > V. > VAIfv6n;z^ c) VAIfven;z > Vz > V.,lo„;z and d) V, < 

ysiow;z- In region a) we keep all seven quantities (V^ is always 
taken from Eq. ISTT l fixed to their analytical values (since all 7 
waves are directed inward), while the number of extrapolated 
variables increases by one each time we cross a critical surface, 
being left with only 4 specified variables in d). In particular P, 
B„, are free to evolve in region d), while only P, in region 
c) and finally only B^ in b). Note that, since the entropy wave is 
always directed inside, at least four out of seven MHD waves are 
associated with a physical boundary condition. 

As far as the rest of the boundaries are concerned, we im- 
pose axisymmetry at the axis and outflow conditions at the top 
z boundary. Moreover, at the outer radial boundary, we apply 
outflow conditions for the ADO, whereas we keep the deriva- 
tive of B^ constant for the simulations of the ASO solution. This 
is particularly important for the latter case where the jet shows 
a high degree of collimation. As a result, a free condition for 
the toroidal component of the magnetic field would cancel the 
poloidal current along the boundary. Hence, the Lorentz force 
would be zero and the uncompensated hoop stress would create 
artificial collimation (as discussed in Zanni et al. 120071 1 . 

On the other hand, when spherical coordinates are adopted, 
axisymmetry holds at the axis and free conditions are imposed 
at the outer radial boundary. At the inner radial boundary and 
at Omax, the correct number of conditions are specified, as pre- 
viously mentioned. Of course, this is achieved according to the 
perpendicular velocities which now are Vr and Ve, respectively. 



4. Results 

4.1. The Analytical Disk Outflow (ADO) solution 
4.1.1. Boundary conditions 

We begin by presenting the results obtained by the correct spec- 
ification of the boundary conditions, as discussed in section ^3.31 
(models IDR, 2DB, 3DB). 

In all cases the initial configurations are maintained and 
are almost identical: density contours and magnetic fieldlines 
are being plotted for model 3DB (same as model 2DB but in 
higher resolution) in Fig. [1] The surfaces corresponding to the 
slow magnetosonic and slow magnetosonic separatrix surface, 
the Alfven surface, and the fast magnetosonic surface of the final 
steady state solution are perfectly coinciding with those of the 
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Fig. 1. Logarithmic density contours (black lines) and poloidal 
magnetic fieldlines (white thin lines) are plotted for the final 
steady-state reached for model 3DB. The critical surfaces of the 
initial configuration are given by the plus signs; going coun- 
terclockwise we find the slow magnetosonic and slow magne- 
tosonic separatrix (which practically coincide), the Alfven, the 
fast magnetosonic, and the fast magnetosonic separatrix surface 
(FMSS, close to the axis). For the final state, the same surfaces 
(apart from the FMSS) are indicated with the thick white lines. 
However, since the analytical equilibrium is preserved in that re- 
gion, they can hardly be distinguished from those predicted by 
the steady-state analysis. On the other hand, the FMSS of the fi- 
nal state is not conical any more but diverges from the axis (its 
position is at the break of the magnetic fieldlines and at the dis- 
continuity of the density). The dashed poloidal fieldline has been 
selected to compute the integrals of motion shown in Fig.|2l 



initial analytical solution. This correspondence was also found 
in GVT06, but it was not so tight: this may be most likely as- 
cribed to the different numerical codes used for the simulations. 
Conversely, the numerically reshaped fast magnetosonic separa- 
trix (FMSS) diverges from the analytical conical position, and, 
as will be discussed later, corresponds to a weak shock which 
can be seen at the density jump and in the break of the poloidal 
fieldlines. 

In Fig. |2l the integrals of motion (Eqs. Q-lO) are being 
shown for models IDR and 2DB respectively, normalized to 
their outer boundary values. The fieldline along which they are 
computed, begins in a region where V, < VAifven;z while it 
crosses the Alfvenic and fast classical critical surfaces. This is 
particularly important since a fieldline with its rooting point al- 
ready in the super- Alfvenic or the super-fast region would give 
much more similar and hence misleading results for the two 
cases. It is clear that the integrals of the final solution deviate 
by less than ~ 1 % when compared to the theoretically expected 
one, with the exception of the entropy integral Q, which shows 
a rather more sensitive behavior. Note an improvement on the 
constancy of those integrals as compared to GVT06, where they 
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Fig. 2. The integrals of motion (Eqs. ||5l-|l8l), normalized to 
unity, are plotted along the poloidal magnetic fieldline (dashed 
line) indicated in Fig. [T] for the models IDR (left) and 2DB 
(right). 
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Fig. 3. The density (left) and total pressure (thermal plus mag- 
netic, right) jumps are being plotted for model 6DH along the 
perpendicular to the shock direction centered at the point x - 
5,y = 33 (indicated with the diamond in Fig.|4|l. 



were found within < 15% of their analytical values, again with 
the exception of Q. 

Therefore, Fig. |2]- along with the fact that both models IDR 
and 2DB give identical final outcomes - suggests that the re- 
sults obtained are not really sensitive to the number of bound- 
ary conditions imposed. The reason for this, is the fact that the 
solution is topologically stable and remains close to the initial 
one. Nevertheless, in order to be physically consistent, the cor- 
rect number of boundary conditions has to be specified always. 

4.1.2. Shock at the FMSS 

We now adopt the setup of 2DB and perform simulations by ef- 
fectively increasing the resolution to examine both the behavior 
and nature of the density and pressure jump observed (models 
4DH, SDH, 6DH). 



Fig. 4. The characteristics (thin solid lines) of the fast magne- 
tosonic waves are plotted in a zoomed region around the shock 
(thick line) for model SDH. The point indicated with a diamond 
is where the quantities of Fig. [3] are plotted. The dashed lines 
are (clockwise) the analytical FMSS, the fast poloidal critical 
surface and the Alfvenic one respectively. 



The discontinuity manifesting in the simulations of the ADO 
models is identified as a weak shock. This conclusion is sup- 
ported by taking into account the negative divergence of the ve- 
locity (thus denoting the compressional nature of the disconti- 
nuity), the jumps appearing on the density and pressure and the 
fact that the gradients steepen with the increase of the resolu- 
tion. The density and the total pressure, thermal plus magnetic 
(P + 12), are plotted across the shock in Fig. [3] Moreover, the 
strength of the shock is found to decrease as we move far away 
from the base, whereas it becomes more and more oblique (see 
Fig. [TJ. Consistent with the shock is the break of the poloidal 
fieldlines, as an effect of the amplification of the tangential com- 
ponent of the poloidal magnetic field. We remark also that the 
jump in the entropy is very small, though, this is expected by the 
almost isothermal conditions assumed for the wind (y - 1.05). 
However, it has been validated (see discussion on solutions with 
different polytropic indexes) that by increasing y — > F the en- 
tropy jump is increasing as well, as expected in adiabatic condi- 
tions. 

In Fig.m selected fast magnetosonic characteristics are plot- 
ted. Note that the cones on the left of the shock do not ever 
cross it but become at best parallel to it. This proves that it cor- 
responds to the numerically readjusted FMSS of the initial exact 
solution. In particular, the model retains its property of a super- 
modified-fast solution, and thus a shock develops to preserve 
the causal disconnection between the flows downstream and up- 
stream. Because the cones constructed by the characteristics of 
the fast waves in the final superfast region never cross the shock 
front, it is not surprising that the sub-modified-fast region is not 
affected at all by the modifications taking place at small polar 
angles. In other words, the new FMSS behaves like a "wall" pre- 
venting the readjustments occuiTing by the extrapolation close 
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Fig. 5. The poloidal currents (loci where tirB^ = const) of 
model SDH are plotted. The arrows indicate the direction of the 
poloidal current density J p. A strong current sheet appears tan- 
gent to the shock. The dashed lines use the notation of Fig.|4] 

to the axis and downstream of the analytical FMSS to affect the 
solution upstream of the shock. Such a result is supported by 
all simulations carried out with the ADO models presented in 
Table □ 

Furthermore, in Fig. |5] we plot the poloidal currents for 
model 5DH. They are counterclockwise upstream and far from 
the FMSS, they change to clockwise very close and upstream 
of the FMSS, and finally they change back to counterclockwise 
downstream of the FMSS. At the FMSS, where the azimuthal 
magnetic field is discontinuous, we have a strong current sheet 
with the surface current heading upwards being tangent to the 
shock. Thus, the currents seem to have a "lightning", or reverse 
"N" shape, with their middle part being parallel to the shock. 
This distribution of the currents, and in particular the direction 
of the resulting J p x force, is consistent with the decolli- 
mation and deceleration that the flow experiences as it passes 
through the shock. Thus, one of the effects of the new FMSS is 
to bend the streamlines away from the z-axis avoiding the over- 
coUimation property of the original analytical solution. The col- 
limation and decoUimation processes that can be derived from 
such a plot are also discussed in GVT06 (see Fig. 6 there). 



Fig. 6. Plot of the initial alfvenic mach number, M, vs log(6') with 
the three extrapolation schemes assumed for models I) 2DB, II) 
IDS and III) 
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4.1 .3. Smoothening the flow near the rotational axis 

We now proceed to study how different types of extending the 
solution up to the axis affect the final outcome of the simulation. 
We adopt the following simple, but diverse enough, extrapola- 
tion schemes for the key functions which are shown in Fig.|6]as 
well: 

{I) (2DB) linear extrap. for < 0„„„ 
II) (IDS) flat extrap. for 6 < 6„i„ (32) 
III) (SDS) smooth extrap. for 0fMss 7.5° 

Notice that case I) is applied to all but IDS and %DS numerical 
models of the ADO solution presented in this paper 



Fig. 7. Logarithmic density contours and magnetic poloidal 
fieldlines for model 8D5 (case III). The symbols are the same 
as in Fig. [T] The results for model IDS (case II) are identical 
with those of models 2DB the morphology of which can be seen 
in Fig.[T] 

Fig. |2] displays the contours of the density for case III) af- 
ter the simulation has reached a steady state, whereas cases I) 
and II) give identical results which can be seen in Fig.[T] These 
first two schemes, which involve modifications only for a small 
polar angle, display very similar features, despite their different 
extrapolation assumption. Case III) is of particular importance 
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Fig. 8. The integrals of motion (Eqs. El -ID) are plotted along m 
at z = 50 for cases I) 2DB (left), II) 7DS (middle) and III) 8DS 
(right). The dashed line represents the initial setup whereas the 
sohd the final steady-state reached. 



since the initial solution is sub-modified-fast, i.e. the whole do- 
main is causally connected. However, the shock is still present 
in the final steady-state reached. Such a result suggests that even 
by starting from an analytical solution that does not cross the 
modified-fast critical point, it will probably self-adapt to an "as- 
trophysically correct" solution. With such a term, we mean that it 
will be consistent with the causal disconnection - of the launch- 
ing and termination regions - expected to exist in the jet phe- 
nomenon. This was somehow anticipated because, for the super- 
modified-fast extrapolation schemes of cases I) and II), the nu- 
merical FMSS is encountered before the analytical one along the 
flow. 

On the other hand, these simulations demonstrate the stabil- 
ity and robustness of the analytical solution since the classical 
poloidal critical surfaces are not readjusted at all for any of these 
cases. We conclude that the final steady-state numerical solu- 
tions do not show any sensitivity neither on the extrapolation 
scheme nor to an initial sub-modified-fast solution as far as the 
criticality condition and the upstream of the shock domain are 
concerned. 

Finally, in order to argue that the extrapolation schemes 
adopted do not coiTespond to any physical inconsistencies, we 
present in Fig.[8]the invariants (Eqs. Q-fS)) along tir at z = 50 
for each case for the initial and final states. Note that case I), 
which is applied in almost all simulations, is particularly inter- 
esting since the inner region is naturally substituted by a lower 
Q mimicking an outflow coming from a slower rotating star. 

4.1 .4. Lower boundary 

We now consider the influence of the choice of z,„,>, to the fi- 
nal steady-state reached since the origin is a point where several 
variables become singular. We construct models 9Dz and IQDz, 



by lowering and increasing z„„„ respectively, and we carry out 
simulations in order to investigate this issue. 

Fig. |9] shows the pressure contours of models 9Dz (left) and 
IQDz (right); the position of the shock is also evident. Note that 
in an iso-density contour plot, the shock would not be clearly 
distinguished because the discontinuity appearing in model 9Dz 
is approximately aligned to the iso-density surfaces, as seen in 
Fig. [T] Even though the classical critical surfaces of both cases 
do not present any deviation from the initial model, the region 
where the shock front develops is considerably different. This 
result can be understood by noting the following. The charac- 
teristics of the fast magnetosonic waves are directly related to 
the formation of the shock. Therefore, the values of the physical 
variables in the extrapolated region control its position. Recall 
that these values are kept fixed since we are in the super-modified 
fast region. So, we speculate that the closer to the origin we set 
the bottom boundary, the stronger these extrapolated values will 
be affected by the steep gradients of the singular origin and the 
more the characteristics defining the shock will deviate from the 
analytical FMSS. 



4.1 .5. Extension to the equator and energetics 

We proceed now to extend GVT06 down to the equator and 
simultaneously examine the effect of different kinds of energy 
input/output. To achieve the former we adopt spherical coordi- 
nates in order to naturally avoid the singularity at the origin, 
where all critical surfaces coincide. The simulations are per- 
formed in the grid with R e [10,90] and 9 e [0,7r/2 - e], 
utilizing R x 9 - [408 x 256] uniform zones, where e ^ 1°. 
We effectively change the way the energy equation is treated by 
implementing and evolving the following 3 setups: i) applying 
y = 1.05 (IIDE), i.e. according to the analytical solution, ii) as- 
suming r - 5/3 (12DG) to examine the solution's behavior to 
an adiabatic evolution and iii) constraining the system under an 
isothermal condition (13DI). 

The left panel of Fig. [10] gives the iso-density surfaces and 
the shape of selected magnetic fieldlines of model 1 IDE. As ex- 
pected, including the equator does not involve any new phenom- 
ena, since the analytical solution is well defined there. However, 
we note that such a computational domain, along with the fixed 
boundary conditions imposed there, is a physically more con- 
sistent choice to describe the MHD outflow, justified by the ge- 
ometric properties of the disk-star system. Finally, notice that 
the position of the shock is consistent with the analysis we per- 
formed with the cylindrical coordinates, with the role played by 
z„„„, being now by /?„„■„. 

On the other hand, interesting results are displayed on the 
middle and right panels of Fig. [10] where selected fieldlines are 
plotted for the initial setup (solid) and the final numerical steady- 
state solution (dashed). The plots indicate that even though the 
polytropic relation is an unavoidable simplification to derive the 
exact solution of VTSTOO, it is found to contribute negligibly on 
the final outcome reached by the numerical evolution. Besides, 
this is not surprising since the necessary distribution of the pres- 
sure needed for constructing such solution is still the same for 
these cases, with only the energetics during time evolution being 
different. Such a pressure is crucial in order to force the Alfvenic 
critical surface to be close to the equator. For a thorough discus- 
sion see Ferreira & Casse ( |2004| l. 
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Fig. 9. Contours of the logarithm of the pressure are being plotted for models 9Dz, i.e. z = 3 (left) and lODz, i.e. z - 12 (right). The 
values increases towards the axis from a minimum value of -4.0 to a maximum of 0.0, while the iso-baric lines differ by 0.1. The 
initial and final poloidal critical surfaces are also indicated with the notation already mentioned in Fig.[T] 




Fig. 10. Logarithmic contours of the density (black lines) are being plotted for models II DE (poly tropic evolution, left), 12DG 
(adiabatic evolution, middle) and 13D/ (isothermal evolution, right). The initial (solid) and final (dashed) fieldlines are sketched as 
well (white lines). 



4.1.6. Long term evolution 

Finally, the structural stability argument has to be made concrete 
by evolving the solution for larger time scales. For that reason 
we constructed model 14Dr, a case identical to 2DB, but with 
its right edge further extended to avoid spurious reflections at the 
boundaries. Moreover, we choose a point in the upstream of the 
shock region, i.e. in the sub-modified-fast domain, and plot the 
deviations of the density from its initial value as a function of 
time. 

It is evident from Fig. [TT] that initially (f < 2) the sub- 
modified-fast solution is being perturbed due to the proper mod- 
ifications at the axis. However, it converges quickly (f < 5) to 
a constant value, roughly ~ 1.5% different from its initial one. 
For the rest of the simulation (5 < f < 220) the steady state is 
perfectly preserved proving its stability. Though, for t > 220, 



boundary effects of the imposed outflow conditions of the right- 
most and upper boundary have propagated throughout the do- 
main and start to artificially affect the solution. It is worth not- 
ing here that the Alfven velocity is of the order of unity at the 
lower central region of the computational box. Hence, the Alfven 
waves have time to propagate many times before the simulation 
is brought to a halt. 

4.2. The Analytical Stellar Outflow (ASO) solution 

The meridionally self-similar solution is in general less com- 
plicated to study, since it does not involve any separatrix in 
the super- Alfvenic region. Therefore, we will mainly investigate 
here its topological stability and its response towards different 
restrictions on the evolution of its energy equation. Notice that 
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Fig. 11. Normalized density deviation from its initial value as a 
function of time at the point {m, z) = (10, 35). 
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Fig, 12. Logarithmic density contours (black lines) are sketched 
for the final time of the simulation for the model ISR which 
does not include a sub-Alfvenic region. The magnetic fieldlines 
are drawn with white lines. The thick dashed one, is chosen to 
plot along the integrals of motion (Eqs. ||5|-|l8l; Fig.[T4li. 



models ISR, 2SL and 3S L evolve in time with the analytically 
derived source term participating throughout the whole compu- 
tational domain. For the rest of the ASO models investigated, the 
details are given in the following. 
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Fig. 13. Selected magnetic fieldlines (white lines) are plotted for 
model 3SL which includes a sub-Alfvenic region. Initial an- 
alytical solution (solid) and final outcome of the simulations 
(dashed). The Alfven critical surface is indicated with a thick 
dotted-dashed line. The gradients of gray represent the logarithm 
of the density. 

4.2.1. Asymptotic configuration 

In Fig. [12] we plot the iso-density surfaces along with selected 
magnetic fieldlines for a super-Alfvenic wind (model ISR). Note 
that in this case we neither have the entropy integral nor the en- 
ergy one, since for the former there is not any polytropic relation 
assumed, while for the latter, an explicit energy source term is 
participating. The initial setup is an exact solution throughout 
the whole computational domain and hence the final time of the 
simulation is arbitrarily chosen to be equivalent to the one se- 
lected for model ISR. The stability of this class of solution can 
be easily seen in Fig. [T4]by the fact that the integrals of motion 
deviate by only < 0.2%. 

The solution stability has also been successfully validated 
throughout all its radial range by adopting a logarithmic grid 
in the radial direction in spherical coordinates (models 2SL 
and 3S L). Note that model 3SL is particularly interesting since 
R € [0.7, 7.0] with Rt - 1, which implies that both sub-Alfvenic 
and super-Alfvenic regions are consistently included inside the 
computational box. Fig. [13] gives the magnetic fieldlines of the 
final numerical solution which remains identical to the initial 
analytical one. Therefore, the application of a logarithmic grid 
allows the approach of the stellar surface, and hence the fixed 
boundaries imposed are physically justified. 

4.2.2. Energetics 

It has already been mentioned that one of the main questions 
posed by the mixing of the two self-similar solutions is on the 
treatment of the energy equation. Although this will be discussed 
extensively in a companion work, we carry out simulations here 
to examine whether the ASO model can reach a steady-state 
when different energy input/output is included. We firstly ad- 
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dress the super-Alfvenic outflow assuming, as in the previous 
ADO solution, y - 1.05 (model \SF) or the adiabatic case 
r = 5/3 (model 55 G). 

The outcome of the simulations is a quasi-steady-state that 
remains very close to the initial solution. With the term quasi- 
steady, we imply the state when the timescale of the evolution of 
the system is much larger as compared to the one of the initial 
more intense readjustments. Especially at the first time steps, 
the solution is strongly perturbed searching for a new equilib- 
rium due to the different energy source terms imposed. Later on, 
after the system relaxes, the outer radial regions are found to 
have remained almost unmodified, whereas, close and along the 
axis the pressure has increased by roughly half an order of mag- 
nitude followed by a slight decrease of the density. This is not 
surprising, since the analytical explicit energy source term, be- 
ing taken into account in model \SR, is of negative sign there, 
i.e. corresponds to energy loss. On the contrary, the polytropic 
case corresponds to a positive input of energy while model 5S G 
to a zero heating/cooling. Hence, the pressure keeps increasing 
in both cases, due to the absence of the needed cooling, un- 
til it reaches a new numerical quasi-equilibrium configuration. 
The quasi-steady-state reached can be judged by Fig. [141 where 
the integrals of motion show deviations of < 3% after ~ 100 
Keplerian rotations at the Alfven radius. This is an expected out- 
come considering that we are in the super-Alfvenic region: we 
know that the efficiency of the thermal driving of the flow is 
concentrated very close to the base, where almost the whole ac- 
celeration occurs. This can be seen in Fig. [15] where we plot 
the energy source term coming from the analytic solution (Eq. 
1301 ). Up to four Alfven radii, the heating distribution decreases 
by ~ 4 orders of magnitude thus proving its crucial role acceler- 
ating the outflow in the inner region. Afterwards, when the flow 
is propagating with its asymptotic speed, the energetics play a 
less important role. 

One of the properties of these solutions are the oscillations 
of the fieldlines. The analytical results predict this possibility 
and we can argue that different energetic processes in the super- 
Alfvenic region amplifies these features. 

The results are totally different if the same type of simula- 
tions, i.e. polytropic or adiabatic evolution, are performed with 
the sub-Alfvenic region included (models 6S L and IS L): in this 
case there is no steady-state reached. On the left of Fig. [16] a 
snapshot of the turbulent evolution is displayed when the poly- 
tropic assumption with y = 1 .05 is applied. In fact, such a model 
is able to drive a sporadic low density outflow around the axis, 
as it can be seen by the velocity vectors. 

Conversely, an adiabatic evolution. Fig. [16] on the right, 
forces the system to collapse towards the star, asymptotically ap- 
proaching a static atmosphere. This is because the ASO model is 
thermally driven and when we impose an adiabatic equation of 
state we effectively switch off' all the heating needed to drive the 
outflow. As expected, the energy processes at the base of merid- 
ionally self-similar winds are crucial for their evolution. 



5. Summary and conclusions 

In this paper we have studied several physical and numerical as- 
pects concerning two classes of the self-similar models, each 
associated with a disk- and a stellar-wind, in the framework 
of the upcoming work that combines them to describe a two- 
component outflow. These analytical solutions (ADO and ASO) 
were appropriately modified, implemented as initial conditions 
and evolved in time. Our main conclusions are the following: 
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Fig. 15. The logarithm of the absolute value of the heating distri- 
bution is being plotted (model 3SL). The thick white line indi- 
cates the surface after which the energy input changes sign and 
becomes cooling. Note that in the inner region, the heating dis- 
tribution gradually decreases by ~ 4 orders of magnitude, while 
in the outer one by only half. 



- The Analytical Disk Outflow ( radially self-similar) solution 
has been successfully validated for its stability and robust- 
ness against several physical and numerical issues. This ar- 
gument holds true even though the analytical solution was in 
many cases significantly modified. We have constructed nu- 
merical models and carried out simulations a) by assuming 
the extreme cases of an isothermal and an adiabatic evolu- 
tion, b) by treating the diverging behavior of the solution at 
the axis with different kinds of extrapolation schemes, mim- 
icking a stellar wind component and c) by changing the size, 
resolution and geometry of the computational box. 
In all cases, the poloidal critical surfaces, with the excep- 
tion of the FMSS, were not readjusted, but rather matched 
perfectly to their initial position. The numerical solution al- 
ways maintained the property of the successful crossing of 
all three critical surfaces producing an outflow causally dis- 
connected from the base. This is achieved with the forma- 
tion of a shock, corresponding to the numerically readjusted 
FMSS. In particular, this shock acts as a "wall" protecting the 
sub-modified-fast magnetosonic regions (source regions of 
the disk wind) from any perturbations taking place due to the 
modification of the models close to the axis (e.g. an effective 
stellar wind). However, the numerically readjusted FMSS 
(shock) does not coincide with the analytical one, with this 
departure being dictated by the respective numerical mod- 
ifications of the models under consideration. A highly sig- 
nificant result is the fact that such a conclusion holds true 
even if we initialize the simulation with a sub-modified-fast 
solution, i.e. a solution with its whole domain causally con- 
nected. We found that, during the simulation, such a numer- 
ical model self-adapts to produce a shock (corresponding to 
the FMSS), hence no information of the downstream region 
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Fig. 16. Logarithmic density contours are being plotted for a polytropic evolution {6SL left) and an adiabatic one (75 L right). The 
velocity vectors overplotted are in the range [4 ■ 10^^,4] and [7 ■ 10"^, 5 ■ 10"'] for the polytropic and adiabatic case, respectively. 
These models never reach a steady-state. 



can travel back to affect the launching region. This implies 
that even MHD outflow solutions, that do not successfully 
cross all three critical points, will probably converge to "as- 
trophysically correct" solutions once evolved in time (see 
also Ferreira 1997 ). 

On the other hand, the study of GVT06 was successfully ex- 
tended down to the equator with the help of simulations us- 
ing spherical coordinates. Furthermore, by adopting different 
assumptions for the energy source terms, it was shown that 
the solution is only slightly and accordingly self-modified 
maintaining all its well defined properties. This is in agree- 
ment with the fact that the ADO solution describes essen- 
tially a magneto-centrifugally accelerated outflow. 
- The Analytical Stellar Outflow (meridionally self-similar) 
solution, which was validated in time-dependent simulations 
for the first time, maintained its well-defined equilibrium as 



expected. Such conclusion is supported by simulations per- 
formed with both the super- and sub-Alfvenic regions in- 
cluded. Quite critical are, contrary to disk winds, the effects 
of the energetics in such thermally driven models. Although 
different assumptions of the energy equation in the super- 
Alfvenic domain did not yield any significant modification 
of the analytical solution, strong variations of the struc- 
ture of the axial outflows are found if modifications of the 
heating/cooling mechanisms occur in the initial accelerat- 
ing region. In particular a polytropic assumption, mimicking 
isothermal conditions, would produce a turbulent weak out- 
flow, while an adiabatic evolution asymptotically reaches a 
static atmosphere. We are tempted to relate the heating inter- 
mittency and even a switching off in such an ASO solution 
with the observed variability of accretion-driven YSO out- 
flows. 
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- All previous statements hold true while being in perfect 
agreement with physically consistent requirements, such as 
specifying the correct type of boundary conditions: a) ac- 
cording to the propagation direction of the MHD waves, b) 
the axisymmetry holding around the axis and c) the con- 
stancy of certain physical variables at both a conical surface 
close to the equatorial plane and a radial one close to the 
origin, implying the presence of an underlying disk and the 
stellar surface, respectively. Therefore, the results can safely 
be trusted since they are not subject to any artificial forcing. 

- Last, but certainly not least, is the fact that almost all mod- 
els of both classes reached a steady- or quasi- steady- state. In 
this context, the upcoming mixing of the two complemen- 
tary classes of self-similar solutions in order to study a two- 
component jet is well founded and promising. The final nu- 
merical model will incorporate both proposed scenarios of a 
pressure-driven outflow (ASO) surrounded by an extended, 
magneto-centrifugally driven disk wind (ADO). This task is 
undertaken in the second paper of this series. 
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